
clear all
cd "______"
log using "HF_Panel_Code.log", replace
*dataset.dta is the wide version of the original dataset (confidential)
use  "dataset.dta", clear

xtset premise_id sasdate

gen month=month(sasdate)
gen year=year(sasdate)
gen dow=dow(sasdate)
gen week=week(sasdate)
gen yearmo=100*year+month
gen yearwk=100*year+week

*keep May-Sep data only
keep if month>4&month<10

forvalues i=1/24{
	gen m_`i'=missing(kwh_`i')
	by premise_id: egen premise_m`i'=sum(m_`i')
}
drop m_*
egen premise_miss=rowtotal(premise_m*)
* premises with less than 1% hourly observations dropped
drop if premise_miss >73
drop premise_m*

egen kwh=rowtotal(kwh_*)
egen temp=rowtotal(temp_*)
egen kwhm=rowmean(kwh_*)
egen tempm=rowmean(temp_*)
drop if kwhm==.

by premise_id: egen numkwh=count(kwh)
by premise_id: egen numtemp=count(temp)
drop if numkwh<300


****Table 1 produced using HF_Panel_Code_long.do

*** Table 2: Summary Statistics
sum kwh kwhm kwh_* tempm temp_*



*** Table B1: Hourly Regressions (Figure 1 below)
mat betahath=J(24,24,.) 
gen temp1=temp_1
forvalues i=1/6{
  local j=`i'
  replace temp1=temp_`j'
  eststo: reghdfe kwh_`i' temp1 , absorb(premise_id yearmo) cluster(yearwk)

  mat betahat=e(b)
  mat betahath[`i',`i']=betahat[1,1]
  reghdfe temp1, absorb(premise_id yearmo) resid
  predict twfe_temp_`i', res
}
esttab using TableB1_1-6.tex, replace b(4)
eststo clear

forvalues i=7/12{
  local j=`i'
  replace temp1=temp_`j'
  eststo: reghdfe kwh_`i' temp1 , absorb(premise_id yearmo) cluster(yearwk)
  mat betahat=e(b)
  mat betahath[`i',`i']=betahat[1,1]
  reghdfe temp1, absorb(premise_id yearmo) resid
  predict twfe_temp_`i', res
}
esttab using TableB1_7-12.tex, replace b(4)
eststo clear

forvalues i=13/18{
  local j=`i'
  replace temp1=temp_`j'
  eststo: reghdfe kwh_`i' temp1 , absorb(premise_id yearmo) cluster(yearwk)
  mat betahat=e(b)
  mat betahath[`i',`i']=betahat[1,1]
  reghdfe temp1, absorb(premise_id yearmo) resid
  predict twfe_temp_`i', res
}
esttab using TableB1_13-18.tex, replace b(4)
eststo clear

forvalues i=19/24{
  local j=`i'
  replace temp1=temp_`j'
  eststo: reghdfe kwh_`i' temp1 , absorb(premise_id yearmo) cluster(yearwk)
  estimates store results_`i'
  mat betahat=e(b)
  mat betahath[`i',`i']=betahat[1,1]
  reghdfe temp1, absorb(premise_id yearmo) resid
  predict twfe_temp_`i', res
}
esttab using TableB1_19-24.tex, replace b(4)
eststo clear
drop temp1

*****Figure 1 
forvalues i=1/24{
gen T`i'=temp_`i'
reghdfe kwh_`i' T`i' , absorb(premise_id yearmo) cluster(yearwk)
estimates store hour`i'
}
*yline at 0.032, which is the average hourly betahat
coefplot hour*, vertical drop(_cons) pstyle(p1) yline(0.032) legend(off) coeflabels(T1 = "1" T2 = "2" T3 = "3" T4="4" T5="5"  T6="6"  T7="7"  T8="8"  T9="9" T10="10" T11="11" T12="12" T13="13" T14="14" T15="15" T16="16" T17="17" T18="18" T19="19" T20="20" T21="21" T22="22" T23="23" T24="24")
graph save Graph "Figure1.gph", replace



/* Table 3, Figures 2 and B1 are produced using the Rcode "Table3_Figure2_FigureB1.R"
 using the output of the csv files "twfecovmat_clean.csv" 
"hourly_betahat_clean.csv" generated using the following code
*/
correlate twfe_temp_*, covariance
return list
matrix twfecovmat=r(C)
putexcel set "twfecovmat_clean.csv", replace
putexcel A2=matrix(twfecovmat)

forvalues i=1/24{
forvalues j=1/24{
if `i'!= `j'{
reghdfe kwh_`i' temp_`j' , absorb(premise_id yearmo) cluster(yearwk)
mat betahat=e(b)
mat betahath[`i',`j']=betahat[1,1]
}
}
}
putexcel set "hourly_betahat_clean.csv", replace
putexcel B2=matrix(betahath)

** computing bethatLF using Eq(30) and the average hourly coefficient betahatbar
reghdfe tempm, absorb(premise_id yearmo) resid
predict twfe_adtemp, res
sum twfe_adtemp
return list
scalar vartwfe_adtemp=r(Var)
scalar list vartwfe_adtemp
matrix C=J(1,24,1)
matrix twfevarvec=vecdiag(twfecovmat)
*betahatLF using Eq(30)
matrix betahatLF=(twfevarvec*(C*betahath)')/(vartwfe_adtemp*(24^2))
matrix list betahatLF
*betahatbar
matrix betahatbar=(C*vecdiag(betahath)')/24
matrix list betahatbar



*** Table B2 (Figure 3A and 3B below)
* Note: generate tempdemean only to use replace in loops

gen tempdemean=tempm
forvalues i=1/6{
  local j=`i'
  gen tempdemean_`i'=temp_`i' - tempm
  replace tempdemean=tempdemean_`i'
  eststo: reghdfe kwh_`i' tempm tempdemean, absorb(premise_id yearmo) cluster(yearwk)

}
esttab using TableB2_01_06.tex, replace b(4)
eststo clear

forvalues i=7/12{
  local j=`i'
  gen tempdemean_`i'=temp_`i' - tempm
  replace tempdemean=tempdemean_`i'
  eststo: reghdfe kwh_`i' tempm tempdemean, absorb(premise_id yearmo) cluster(yearwk)
}
esttab using TableB2_07_12.tex, replace b(4)
eststo clear

forvalues i=13/18{
  local j=`i'
  gen tempdemean_`i'=temp_`i' - tempm
  replace tempdemean=tempdemean_`i'
  eststo: reghdfe kwh_`i' tempm tempdemean, absorb(premise_id yearmo) cluster(yearwk)
}
esttab using TableB2_13_18.tex, replace b(4)
eststo clear

forvalues i=19/24{
  local j=`i'
  gen tempdemean_`i'=temp_`i' - tempm
  replace tempdemean=tempdemean_`i'
  eststo: reghdfe kwh_`i' tempm tempdemean, absorb(premise_id yearmo) cluster(yearwk)
}
esttab using TableB2_19_24.tex, replace b(4)
drop T*

forvalues i=1/24{
  gen M`i'=tempm
  gen T`i'=temp_`i'-tempm
  reghdfe kwh_`i' M`i' T`i', absorb(premise_id yearmo) cluster(yearwk)
  estimates store hourcm`i'
}

** Figure 3A &B
*Figure 3A
coefplot hourcm*, vertical drop(_cons T*)  pstyle(p1) legend(off) coeflabels(M1 = "1" M2 = "2" M3 = "3" M4="4" M5="5"  M6="6"  M7="7"  M8="8"  M9="9" M10="10" M11="11" M12="12" M13="13" M14="14" M15="15" M16="16" M17="17" M18="18" M19="19" M20="20" M21="21" M22="22" M23="23" M24="24")
graph save Graph "Figure3A.gph", replace
*Figure 3B
coefplot hourcm*, vertical drop(_cons M*)  pstyle(p1) legend(off) coeflabels(T1 = "1" T2 = "2" T3 = "3" T4="4" T5="5"  T6="6"  T7="7"  T8="8"  T9="9" T10="10" T11="11" T12="12" T13="13" T14="14" T15="15" T16="16" T17="17" T18="18" T19="19" T20="20" T21="21" T22="22" T23="23" T24="24")
graph save Graph "Figure3B.gph", replace

*Table 5 in paper/Tables B3-6 in online appendix
eststo clear
forvalues i=1/24{
  eststo: reghdfe kwh_`i' temp_`i', absorb(premise_id yearmo) cluster(yearwk)
  eststo: reghdfe kwh_`i' temp_`i', absorb(i.premise_id#i.year yearmo) cluster(yearwk)

  eststo: reghdfe kwh_`i' temp_`i', absorb(premise_id i.yearmo#i.dow) cluster(yearwk)
  eststo: reghdfe kwh_`i' temp_`i', absorb(premise_id yearmo dow) cluster(yearwk)
  eststo: reghdfe kwh_`i' temp_`i', absorb(premise_id year month dow) cluster(yearwk)
  esttab using Table5_`i'.tex, replace b(4)
  eststo clear
  }

* Table B7: FD estimator
gen kwhmdiff=kwhm-L.kwhm
gen tempmdiff=tempm-L.tempm
eststo clear
eststo: reghdfe kwhmdiff tempmdiff, noabsorb cluster(yearwk)
eststo: reghdfe kwhmdiff tempmdiff, absorb(yearmo) cluster(yearwk)
esttab using TableB7.tex




log close
